
	SUBROUTINE Shoreln3(IMAX,JSHORE,YOFF,H,Htop,Hbot,DELTAY,DeltaH,
	. tanALPHA,TA,Hav,HAA,HBB,HA,HB,Ylamax,Ylbmax,FMa,FMb,Coef,
	. YSHORE,Hold)
c
c	 This subroutine causes the shoreline to advance or retreat in sub-grid
c      increments.  It is based on the Cowell and Bruun idea that because the
c      time-scale of surf zone and uppermost shoreface profiles is very fast
c      they can be represented as having a fixed shape in long-term simulations
c      with the understanding that they are surrounded by an envelope of profiles
c      that represent the day-to-day conditions.  
c
c      The height change from subroutine SurfZone is redistributed by this
c      subroutine between the beach/surf zone cell (JSHORE(I)) and the 
c      adjoining cells in the I-direction because the fixed profile shape 
c	 (represented by a sloped straight line) slide in, or out (not up
c      or down) as the JSHORE(I) cell empties or fills. As this profile
c      'slides' in or out sand volume is subtracted or added to the
c      adjacent cells accordingly. The cell seaward of JSHORE(I) also
c      adjusts coming up as the beach/surf zone profile advances or
c      going down in the retreating condition.
c
c	 When this 'sliding' of the profile causes the shoreline to move 
c      one of the adjacent cells the values of JSHORE(I) is changed
c      accordingly.
c
	Dimension H(100,100),Hupld(100),JSHORE(100),Yoff(100),
	. SXser(1,100),SYser(1,100),H(150,150),YSHORE(100),
	. Hold(100)
c
c
	DO 800 I=1,IMAX
c
	Hinc = (H(I,JSHORE(I)))-Hold(I)
c
c
c
c	Adjust the upper shoreface cell (i.e. JSHORE(I)-1) so that it comes up
c	 Hbot as the shoreline advances or drops to the adjacent cell heights as
c      the shoreline retreats.  The concept here is that the response time of 
c      surf zone is so fast that it can be thought of as a fixed profile 
c      surrounded by an envelope of maximum and minimum shapes. Next cell
c	 seaward is the upper shoreface which tracks the surf zone in a damped
c	 manner.
c
c	 The movement of this upper shoreface cell depends first on whether the
c      surf zone is advancing or retreating landward.
c
	If(Hinc.GE.0.)Goto 14 ! advancing case transfers
c	  eroding case falls through

	Htarget = (Hbot + (H(I,(JSHORE(I)-2))))/2.
	Halter = Coef*(Htarget-(H(I,(JSHORE(I)-1))))
	IF(Halter.GT.Hinc)GOTO12 !If Halter is more (neg) than Hinc,limit it
	Halter = Hinc
c
   12	Hinc = Hinc - Halter
	H(I,(JSHORE(I)-1)) = H(I,(JSHORE(I)-1)) + Halter
	GOTO18
c     The following is the advancing case
   14	Halter = Coef*(Hbot - H(I,(JSHORE(I)-1)))
	IF(Halter.LT.Hinc)GOTO16 !If Halter is more than Hinc,limit it
	Hchange = Hinc
  16	Hinc = Hinc - Halter
	H(I,(JSHORE(I)-1)) = H(I,(JSHORE(I)-1)) + Halter
c
   18 Continue !If(Hinc.EQ.0.)GOTO800
c
c     Now that the upper shoreface cell is adjusted the surf zone changes
c	are calculated. 
c
c     When converting to the subroutine we will need to compute a Hinc from
c     the difference of Hold(I) and H(I,JSHORE(I)). At this time we set this 
c     difference as an input and make the following statement
c
c	 
	 H(I,JSHORE(I)) = Hold(I) + Hinc
c
c	
c    Test whether the shoreline is advancing seaward or retreating landward
c
	IF(H(I,(JSHORE(I))).LT.Hold(I))GOTO24
c
c        The advancing case falls through
c
	IF((H(I,(JSHORE(I))).GT.Hav).AND.((Hold(I)).GE.Hav))GOTO23
	IF((H(I,(JSHORE(I))).GE.Hav).AND.((Hold(I)).LE.Hav))GOTO22
c
c         The advancing <Hav,<Hav case falls through
c
	Ylbold = Ylbmax*((Hold(I)-Hav)/(HA-Hav))
	Ylbnew = Ylbmax*((H(I,(JSHORE(I)))-Hav)/(HA-Hav))
c
	VolBold = (0.5)*((Ylbold**2)*tanAlPHA)	   
	VolBnew = (0.5)*((Ylbnew**2)*tanALPHA)	   
c
	VolDiff = VolBold - VolBnew
	hDiff1 = FMb*((Ylbold/Ylbmax))*(VolDiff/DeltaY)	  !!!old for new & max)**2
	hDiff2 = FMb*((VolDiff/DeltaY)-hDiff1)
c
	H(I,(JSHORE(I))) = H(I,(JSHORE(I))) - hDiff2
	H(I,(JSHORE(I)+1)) = H(I,(JSHORE(I)+1)) + hDiff1
	GOTO40
c
   22	Continue ! This is the advancing Straddle Hav case
c
	Ylbold = -(Ylbmax*((Hold(I)-Hav)/(Hav-HA)))
	Ylanew = - (Ylamax*((H(I,(JSHORE(I)))-Hav)/(HB-Hav)))
c
	VolBold = (0.5)*((Ylbold**2)*tanAlPHA)
	VolAnew = (0.5)*((Ylanew**2)*tanALPHA)
c
	H(I,(JSHORE(I)+1)) = H(I,(JSHORE(I)+1))+(VolBold/DeltaY)
	H(I,(JSHORE(I))) = H(I,(JSHORE(I))) - ((VolBold+VolAnew)/DeltaY)
	H(I,(JSHORE(I)-1)) = H(I,(JSHORE(I)-1))	+ (VolAnew/DeltaY)
	GOTO40
c
   23	 Continue ! This is the advancing >Hav,>Hav case
c
	Ylaold = - (Ylamax*((Hold(I)-Hav)/(HB-Hav)))
 	Ylanew = - (Ylamax*((H(I,(JSHORE(I))))-Hav)/(HB-Hav))
c
	VolAold = (0.5)*((Ylaold**2)*tanALPHA)
 	VolAnew = (0.5)*((Ylanew**2)*tanALPHA)
c
	diffVol = VolAnew-VolAold
	hDiff1 = ABS(FMa*((Ylanew/Ylamax))*(diffVol/DeltaY))		
	hDiff2 = ABS(FMa*((VolDiff/DeltaY)-hDiff1))
c
	H(I,(JSHORE(I))) = H(I,(JSHORE(I))) - hDiff2
	H(I,(JSHORE(I)-1)) = H(I,(JSHORE(I)-1))	+ hDiff1
	GOTO40
c
c
c      There follows the retreating cases
c
   24	IF((H(I,(JSHORE(I))).GT.Hav).AND.((Hold(I)).GE.Hav))GOTO27
 	IF((H(I,(JSHORE(I))).LE.Hav).AND.((Hold(I)).GE.Hav))GOTO26
c
c	 The retreating <Hav,<Hav case falls through
c
	Ylbold = Ylbmax*((Hold(I)-Hav)/(HA-Hav))
 	Ylbnew = Ylbmax*((H(I,(JSHORE(I)))-Hav)/(HA-Hav))
c
	VolBold = (0.5)*((Ylbold**2)*tanAlPHA)
	VolBnew = (0.5)*((Ylbnew**2)*tanALPHA)
c
	VolDiff =  VolBnew - VolBold 
	hDiff1 = FMb*((Ylbold/Ylbmax))*(VolDiff/DeltaY)	!!!old for new	!!!here
	hDiff2 = FMb*((VolDiff/DeltaY) - hDiff1)
c
	H(I,(JSHORE(I))) = H(I,(JSHORE(I))) + hDiff2
	H(I,(JSHORE(I)+1)) = H(I,(JSHORE(I)+1)) - hDiff1
	GOTO40
c
c
   26	Continue ! This is the retreating Straddle Hav case
c
	Ylaold =(Ylbmax*((Hold(I)-Hav)/(HB-Hav)))
	Ylbnew =(Ylamax*(H(I,(JSHORE(I)))-Hav)/(HA-Hav))
c
	VolAold = (0.5)*((Ylaold**2)*tanAlPHA)
	VolBnew = (0.5)*((Ylbnew**2)*tanALPHA)
c
	H(I,(JSHORE(I)+1)) = H(I,(JSHORE(I)+1))- (VolBnew/DeltaY)
	H(I,(JSHORE(I))) = H(I,(JSHORE(I))) + ((VolAold+VolBnew)/DeltaY)
	H(I,(JSHORE(I)-1)) = H(I,(JSHORE(I)-1))	- (VolAold/DeltaY)
	GOTO40
c
   27	 Continue ! This is the retreating >Hav,>Hav case
c
	Ylaold = - (Ylamax*((Hold(I)-Hav)/(HB-Hav)))
 	Ylanew = - (Ylamax*((H(I,(JSHORE(I))))-Hav)/(HB-Hav))
c
	VolAold = (0.5)*((Ylaold**2)*tanALPHA)
 	VolAnew = (0.5)*((Ylanew**2)*tanALPHA)
c
	diffVol = VolAold - VolAnew
	hDiff1 = ABS(FMa*((Ylaold/Ylamax))*(diffVol/DeltaY))	!!!here 
	hDiff2 = ABS(FMa*((VolDiff/DeltaY)-hDiff1))
c
	H(I,(JSHORE(I))) = H(I,(JSHORE(I))) + hDiff2
	H(I,(JSHORE(I)-1)) = H(I,(JSHORE(I)-1))	- hDiff1
	GOTO40
c
   40	Continue
c
c	Trap the overfilling condition
c
	IF(H(I,(JSHORE(I))).LE.HB)GOTO42
	Delta = H(I,(JSHORE(I)))-HB
	H(I,(JSHORE(I)-1)) = H(I,(JSHORE(I)-1)) + Delta
	H(I,(JSHORE(I))) = HB
c	
c     The shoreline is tracked with the cell-referenced value YPshore which is
c	measured in the plus-Y direction from the right (i.e.offshore) cell face  
c
   42	YPShore = DeltaY*(1.-((H(I,(JSHORE(I)))-HA)/(HB-HA)))
c
c	 The shoreline is also tracked in real grid space by YOFF(I)
c
	YOFF(I) = ((JSHORE(I)*DeltaY) - (DeltaY/2.)) + YPShore
c
c      Determine if a cell shift is needed 
c
	IF(YPShore.GT.0.)GOTO45
	IF(H(I,(JSHORE(I)-1)).GE.HA)GOTO44
	Delta = HA -  H(I,(JSHORE(I)-1))
	H(I,(JSHORE(I))) = H(I,(JSHORE(I))) - Delta
	H(I,(JSHORE(I)-1)) = HA
	GOTO50
   44	JSHORE(I) = JSHORE(I) - 1
	GOTO50
   45	IF(YPShore.LE.DeltaY)GOTO50
	JSHORE(I) = JSHORE(I) + 1
c
c    Store the JSHORE(I) depth for the next time step
c
   50	Hold(I) = H(I,JSHORE(I))
c
c
  800 CONTINUE
c
c
	RETURN
	END